If you are new yo R, make sure to read these info boxes.
First load some dependencies and create the output directory:
In R, packages can be loaded using library(). If a package is not installed, you can install it from CRAN using install.packages(). The dplyr package is a commonly used package for data manipulation, readxl is required for reading the Excel file.
list.files("./dataset", full.names = "TRUE")## [1] "./dataset/metadata.txt" "./dataset/samples.xlsx"
## [3] "./dataset/seqtab.txt" "./dataset/sequences.fasta"
## [5] "./dataset/taxonomy.txt"In a relative file path, .. indicates the parent directory.
../dataset/seqtab.txt contains the ASV table, so it has one row per ASV, and the number of reads in a sample in different columns.
seqtab <- read.table("./dataset/seqtab.txt", sep = "\t", header = TRUE)
rmarkdown::paged_table(seqtab)read.table() reads a delimited text file to a data frame. sep = “ means that the file is tab delimited.
../dataset/taxonomy.txt contains a taxon name for each ASV.
taxonomy <- read.table("./dataset/taxonomy.txt", sep = "\t", header = TRUE)
paged_table(taxonomy)These names originate from the reference database and will have to be matched to WoRMS later.
We also have an Excel file with sample info.
samples <- read_excel("./dataset/samples.xlsx")
samples## # A tibble: 2 × 11
## name size event_begin area_name area_longitude area_latitude
## <chr> <dbl> <chr> <chr> <dbl> <dbl>
## 1 EE0493 1450 24/04/2023 Ile esprit 46.2 9.43
## 2 EE0495 1500 02/04/2023 Settlement beach 46.2 9.40
## # ℹ 5 more variables: area_uncertainty <dbl>, parent_area_name <chr>,
## # dna <dbl>, depth <dbl>, temperature <dbl>At this point we could start quality control on the individual tables, but if we first join and map the tables to Darwin Core occurrence terms, the quality control code will be easier to read.
Let’s start with the sample table. This table has sample identifiers, time, coordinates, coordinate uncertainty, locality, and higher geography which can all be mapped to Darwin Core. Keep the dna field for later.
event <- samples %>%
select(
eventID = name,
materialSampleID = name,
eventDate = event_begin,
locality = area_name,
decimalLongitude = area_longitude,
decimalLatitude = area_latitude,
coordinateUncertaintyInMeters = area_uncertainty,
higherGeography = parent_area_name,
minimumDepthInMeters = depth,
maximumDepthInMeters = depth,
sampleSizeValue = size,
dna,
temperature
) %>%
mutate(sampleSizeUnit = "ml")
event## # A tibble: 2 × 14
## eventID materialSampleID eventDate locality decimalLongitude decimalLatitude
## <chr> <chr> <chr> <chr> <dbl> <dbl>
## 1 EE0493 EE0493 24/04/2023 Ile espr… 46.2 9.43
## 2 EE0495 EE0495 02/04/2023 Settleme… 46.2 9.40
## # ℹ 8 more variables: coordinateUncertaintyInMeters <dbl>,
## # higherGeography <chr>, minimumDepthInMeters <dbl>,
## # maximumDepthInMeters <dbl>, sampleSizeValue <dbl>, dna <dbl>,
## # temperature <dbl>, sampleSizeUnit <chr>The code above uses several dlyr functions. select() selects and optionally renames a set of columns from the dataframe. mutate() creates a new column. %>% is the pipe operator which is used to string functions together.
Next is the ASV table. This table is in a wide format with ASVs as rows and samples as columns. We will convert this to a long format, with one row per occurrence and the number of sequence reads as organismQuantity. We will use the sample identifier as eventID and the combination of sample identifier and ASV number as the occurrenceID.
To do from a wide to a long table, use the gather() function from the tidyr package. paste0() is used to combine character strings.
#library(tidyr)
occurrence <- seqtab %>%
gather(eventID, organismQuantity, 2:3) %>%
filter(organismQuantity > 0) %>%
mutate(
occurrenceID = paste0(eventID, "_", asv),
organismQuantityType = "sequence reads"
)
paged_table(occurrence)We can now add the taxonomic names to our occurrence table.
taxonomy <- taxonomy %>%
select(asv, verbatimIdentification = taxonomy)occurrence <- occurrence %>%
left_join(taxonomy, by = "asv")
paged_table(occurrence)left_join() joins two dataframes by matching columns. The by argument specifies the columns to match on.
occurrence <- event %>%
left_join(occurrence, by = "eventID")
paged_table(occurrence)Populate samplingProtocol with a link the the eDNA Expeditions protocol.
occurrence$samplingProtocol <- "https://github.com/BeBOP-OBON/UNESCO_protocol_collection"Let’s first match the taxa with WoRMS. This can be done using the obistools package. Before matching with WoRMS we will remove underscores from the scientific names.
taxon_names <- stringr::str_replace(occurrence$verbatimIdentification, "_", " ")Now match the names, this can take a few minutes.
matched <- obistools::match_taxa(taxon_names, ask = FALSE) %>%
select(scientificName, scientificNameID)## 433 names, 0 without matches, 11 with multiple matches#The solution can be found in this file:
#matched <- read.table("./solutions/matched_results.txt", sep = "\t", header = T)
paged_table(matched)occurrence <- bind_cols(occurrence, matched)
paged_table(occurrence)non_matches <- occurrence %>%
filter(is.na(scientificNameID)) %>%
group_by(verbatimIdentification) %>%
summarize(n = n()) %>%
arrange(desc(n))
write.table(non_matches, file = file.path(output_dir, "nonmatches.txt"), sep = "\t", row.names = FALSE, na = "", quote = FALSE)
paged_table(non_matches)Normally we have to resolve these names one by one, but for this exercise we will just fix the most common errors. For example, records annotated as eukaryotes can be populated with scientificName Incertae sedis and scientificNameID urn:lsid:marinespecies.org:taxname:12.
occurrence <- occurrence %>%
mutate(
scientificName = case_when(verbatimIdentification %in% c("Eukaryota", "undef_Eukaryota", "") ~ "Incertae sedis", .default = scientificName),
scientificNameID = case_when(verbatimIdentification %in% c("Eukaryota", "undef_Eukaryota", "") ~ "urn:lsid:marinespecies.org:taxname:12", .default = scientificNameID)
)OBIS will automatically link higher taxonomic levels based on the worms IDS. However, with the following workflow you can add them also.
#library(dplyr)
#library(purrr)
dummy_data <- occurrence %>%
select(scientificName, scientificNameID) %>%
mutate(aphiaid = as.numeric(stringr::str_extract(scientificNameID, "\\d+$")))
taxonomy_worrms <- map(unique(dummy_data$aphiaid[!is.na(dummy_data$aphiaid)]), worrms::wm_record) %>%
bind_rows() %>%
select(AphiaID, kingdom, phylum, class, order, family, genus, scientificname, rank)
dummy_data <- dummy_data %>%
left_join(taxonomy_worrms, by = c("aphiaid" = "AphiaID"))
occurrence <- bind_cols(occurrence, dummy_data %>% select(kingdom, phylum, class, order, family, genus, rank))paged_table(occurrence)